1 Setup

1.1 Load libraries

library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")
Data last updated on 2024-05-03

1.2 Handy functions not in spat1f

trend_by_pre_wt <- function(model, term){
  emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>% 
  pairs()
}

contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
  type <- match.arg(type)
  f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
  emmeans::emmeans(model, 
                 f, 
                 var = term, 
                 at = list(cat_pre_wt_log_scale = at))
}

check_mod <- function(model){
 DHARMa::simulateResiduals(model) %>% 
    plot()
  performance::posterior_predictive_check(model)
}

1.3 Quick cleaning

# Main data.frame for analyses
d <- ref_data %>% 
  filter(error == 0) %>% 
  mutate(
    cat_pre_wt_log = log(cat_pre_wt),
    cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
    cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
    has_var = ifelse(var_trt != "constant", 1, 0),
    mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
    mean_toxic_scale = as.numeric(scale(mean_toxic)),
    area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
    var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
    on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic, 
                                                        trans = "emp", 
                                                        nudge.method = "none", 
                                                        na.action = "ignore"))),
    select_scale = as.numeric(scale(less_toxic)),
    ava_qual_scale = as.numeric(scale(ava_qual)),
    ava_qual_logit_scale = as.numeric(scale(adjust_prop(ava_qual, 
                                                        trans = "emp", 
                                                        nudge.method = "none", 
                                                        na.action = "ignore"))),
    sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
    cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
    RGR_scale = as.numeric(scale(RGR)), 
    sl_mean_pred = shape * scale
  )

# Subset of data.frame for SEM
d2 <- d %>% 
  filter(
    var_trt != "constant"
  ) %>% 
  filter(
    !is.na(area_herb_log_scale) & 
      !is.na(var_toxic_12_scale) & 
      !is.na(mean_toxic_conc_scale) &
    !is.na(on_toxic_logit_scale) &
     !is.na(select_scale) &
      !is.na(ava_qual_scale) &
     !is.na(sl_mean_obs) &
      !is.na(RGR)) %>% 
  mutate(
    var_high = ifelse(var_trt == "high_var", 1, 0),
    beta_red = ifelse(as.character(beta) == 5, 1, 0),
    beta_white = ifelse(as.character(beta) == 0, 1, 0),
    beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
    beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
    var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
    beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
  )

2 Herbivore Performance

2.1 RGR

2.1.1 Model summary

rgr_m <- glmmTMB(
  RGR ~ 
    (var_trt + beta) * cat_pre_wt_log_scale +
    I(cat_pre_wt_log_scale^2) +  # residual plot show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(
      var_trt != "constant"
    ) 
); summary(rgr_m)
 Family: gaussian  ( identity )
Formula:          
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
  -879.2   -849.4    450.6   -901.2      100 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev.
 session_id (Intercept) 3.979e-06 0.001995
 Residual               1.603e-05 0.004004
Number of obs: 111, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 1.6e-05 

Conditional model:
                                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          0.0162504  0.0013457  12.075  < 2e-16 ***
var_trtlow_var                      -0.0034872  0.0007740  -4.506 6.62e-06 ***
beta0                               -0.0003065  0.0009664  -0.317  0.75112    
beta5                               -0.0013545  0.0009173  -1.477  0.13978    
cat_pre_wt_log_scale                -0.0048759  0.0011991  -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2)           -0.0021953  0.0006731  -3.262  0.00111 ** 
var_trtlow_var:cat_pre_wt_log_scale  0.0025873  0.0007897   3.276  0.00105 ** 
beta0:cat_pre_wt_log_scale           0.0011924  0.0009491   1.256  0.20897    
beta5:cat_pre_wt_log_scale           0.0024731  0.0009710   2.547  0.01087 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(rgr_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: RGR
                                Chisq Df Pr(>Chisq)    
(Intercept)                  145.8141  1  < 2.2e-16 ***
var_trt                       20.3008  1  6.617e-06 ***
beta                           2.3963  2   0.301758    
cat_pre_wt_log_scale          16.5348  1  4.777e-05 ***
I(cat_pre_wt_log_scale^2)     10.6374  1   0.001108 ** 
var_trt:cat_pre_wt_log_scale  10.7334  1   0.001052 ** 
beta:cat_pre_wt_log_scale      6.4936  2   0.038899 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1.2 Post-hoc contrasts

trend_by_pre_wt(rgr_m, "var_trt")
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00259 0.00079 100  -3.276  0.0014

Results are averaged over the levels of: beta 
contrast_by_pre_wt(rgr_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var  0.014224 0.00355 100  0.007178  0.02127
 low_var   0.005562 0.00312 100 -0.000635  0.01176

cat_pre_wt_log_scale =  2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var -0.000393 0.00208 100 -0.004512  0.00373
 low_var   0.001295 0.00238 100 -0.003421  0.00601

Results are averaged over the levels of: beta 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var  0.00866 0.00181 100   4.790  <.0001

cat_pre_wt_log_scale =  2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00169 0.00171 100  -0.988  0.3256

Results are averaged over the levels of: beta 
trend_by_pre_wt(rgr_m, "beta")
 contrast         estimate       SE  df t.ratio p.value
 (beta-5) - beta0 -0.00119 0.000949 100  -1.256  0.4232
 (beta-5) - beta5 -0.00247 0.000971 100  -2.547  0.0329
 beta0 - beta5    -0.00128 0.000940 100  -1.363  0.3644

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(rgr_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5    0.012890 0.00382 100  0.005313  0.02047
 0     0.010199 0.00333 100  0.003591  0.01681
 5     0.006589 0.00319 100  0.000255  0.01292

cat_pre_wt_log_scale =  2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5   -0.001439 0.00226 100 -0.005922  0.00304
 0     0.000639 0.00231 100 -0.003940  0.00522
 5     0.002153 0.00257 100 -0.002946  0.00725

Results are averaged over the levels of: var_trt 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0  0.00269 0.00221 100   1.218  0.4456
 (beta-5) - beta5  0.00630 0.00223 100   2.820  0.0158
 beta0 - beta5     0.00361 0.00211 100   1.711  0.2062

cat_pre_wt_log_scale =  2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0 -0.00208 0.00205 100  -1.015  0.5689
 (beta-5) - beta5 -0.00359 0.00206 100  -1.746  0.1935
 beta0 - beta5    -0.00151 0.00211 100  -0.719  0.7529

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.1.3 Misc stats

r2(rgr_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.555
     Marginal R2: 0.444

2.1.4 Check model

check_mod(rgr_m)

2.1.5 Plot

g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = expression(beta), fill = expression(beta))

g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) + 
  geom_line(aes(color = var_trt), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("fill","color"), 
                       label = c("High", "Low")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = "Variation", fill = "Variation")


ggarrange(g1, g2 + labs(y = ""))

2.2 Time to pupation

2.2.1 Model summary

pupt_m_full <- glmmTMB(
  time_to_pupation ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m_full)
 Family: Gamma  ( log )
Formula:          time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +  
    I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   338.1    366.3   -158.0    316.1       85 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0254 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          2.180302   0.036853   59.16  < 2e-16 ***
var_trtlow_var                       0.102145   0.033998    3.00  0.00266 ** 
beta0                               -0.001845   0.040823   -0.05  0.96395    
beta5                                0.059543   0.041181    1.45  0.14821    
cat_pre_wt_log_scale                -0.428662   0.032643  -13.13  < 2e-16 ***
I(cat_pre_wt_log_scale^2)           -0.081267   0.019069   -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale  0.026518   0.038876    0.68  0.49516    
beta0:cat_pre_wt_log_scale          -0.030363   0.043195   -0.70  0.48210    
beta5:cat_pre_wt_log_scale          -0.126801   0.047095   -2.69  0.00709 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                                 Chisq Df Pr(>Chisq)    
(Intercept)                  3500.2207  1  < 2.2e-16 ***
var_trt                         9.0264  1   0.002661 ** 
beta                            2.7351  2   0.254736    
cat_pre_wt_log_scale          172.4448  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)      18.1614  1  2.029e-05 ***
var_trt:cat_pre_wt_log_scale    0.4653  1   0.495162    
beta:cat_pre_wt_log_scale       7.5589  2   0.022836 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupt_m <- glmmTMB(
  time_to_pupation ~ 
    var_trt + beta * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m)
 Family: Gamma  ( log )
Formula:          
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   336.5    362.2   -158.3    316.5       86 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0255 

Conditional model:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 2.180555   0.036937   59.03  < 2e-16 ***
var_trtlow_var              0.107447   0.033197    3.24  0.00121 ** 
beta0                      -0.003578   0.040842   -0.09  0.93019    
beta5                       0.059536   0.041271    1.44  0.14915    
cat_pre_wt_log_scale       -0.420147   0.030221  -13.90  < 2e-16 ***
I(cat_pre_wt_log_scale^2)  -0.083365   0.018881   -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011   0.042849   -0.61  0.54383    
beta5:cat_pre_wt_log_scale -0.126834   0.047183   -2.69  0.00719 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                              Chisq Df Pr(>Chisq)    
(Intercept)               3485.0791  1  < 2.2e-16 ***
var_trt                     10.4759  1   0.001209 ** 
beta                         2.8053  2   0.245947    
cat_pre_wt_log_scale       193.2743  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)   19.4946  1  1.009e-05 ***
beta:cat_pre_wt_log_scale    7.6890  2   0.021397 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.2 Post-hoc contrasts

contrast_by_pre_wt(pupt_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  2.808 0.0769 Inf     2.657      2.96
 low_var   2.915 0.0818 Inf     2.755      3.08

cat_pre_wt_log_scale =  2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  0.924 0.0672 Inf     0.792      1.06
 low_var   1.031 0.0725 Inf     0.889      1.17

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

cat_pre_wt_log_scale =  2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
trend_by_pre_wt(pupt_m, "beta")
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0    0.026 0.0428 Inf   0.607  0.8163
 (beta-5) - beta5    0.127 0.0472 Inf   2.688  0.0197
 beta0 - beta5       0.101 0.0473 Inf   2.132  0.0834

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(pupt_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    2.741 0.0941 Inf     2.557      2.93
 0     2.790 0.0971 Inf     2.599      2.98
 5     3.054 0.1088 Inf     2.841      3.27

cat_pre_wt_log_scale =  2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    1.061 0.0845 Inf     0.895      1.23
 0     1.005 0.0805 Inf     0.847      1.16
 5     0.866 0.0917 Inf     0.687      1.05

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0  -0.0484 0.1018 Inf  -0.476  0.8828
 (beta-5) - beta5  -0.3132 0.1133 Inf  -2.763  0.0158
 beta0 - beta5     -0.2648 0.1136 Inf  -2.331  0.0516

cat_pre_wt_log_scale =  2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0   0.0556 0.0875 Inf   0.636  0.8005
 (beta-5) - beta5   0.1941 0.0915 Inf   2.122  0.0854
 beta0 - beta5      0.1385 0.0928 Inf   1.492  0.2945

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.2.3 Misc stats

r2(pupt_m, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.881

2.2.4 Check model

check_mod(pupt_m)

2.2.5 Plot

g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log2") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
       color = expression(beta), fill = expression(beta))

g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>% 
  ggplot(aes(x = var_trt, y = yhat)) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = var_trt, y = time_to_pupation, color = var_trt),
    position = position_jitter(width = 0.2),
    size = 3, alpha = 0.8
  ) + 
  geom_pointrange(
    aes(ymax = upper, ymin = lower), 
    color = "black",
    size = 1, linewidth = 2, shape = 1,
  ) +
  theme_bw(base_size = 15) + 
  scale_y_continuous(trans = "log2") +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("color"), 
                       label = c("High", "Low")) + 
  labs(x = "Variation", y = "Time to pupation (days)",
       color = "Variation")


ggarrange(g3, g4 + labs(y = ""))

2.3 Eclosure

2.3.1 Model summary

eclosure_m_full <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m_full)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   185.3    210.9    -83.6    167.3      119 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1661   0.4076  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)
(Intercept)                          0.53530    0.41050   1.304    0.192
var_trtlow_var                       0.01707    0.37507   0.046    0.964
beta0                               -0.21544    0.46265  -0.466    0.641
beta5                               -0.50749    0.45596  -1.113    0.266
cat_pre_wt_log_scale                 0.10452    0.36089   0.290    0.772
var_trtlow_var:cat_pre_wt_log_scale  0.17526    0.36489   0.480    0.631
beta0:cat_pre_wt_log_scale           0.24778    0.43553   0.569    0.569
beta5:cat_pre_wt_log_scale           0.20159    0.44699   0.451    0.652
car::Anova(eclosure_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: eclosed
                              Chisq Df Pr(>Chisq)
(Intercept)                  1.7004  1     0.1922
var_trt                      0.0021  1     0.9637
beta                         1.2502  2     0.5352
cat_pre_wt_log_scale         0.0839  1     0.7721
var_trt:cat_pre_wt_log_scale 0.2307  1     0.6310
beta:cat_pre_wt_log_scale    0.3659  2     0.8328
eclosure_m <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   179.9    197.0    -84.0    167.9      122 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1676   0.4094  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)           0.551286   0.413578   1.333    0.183
var_trtlow_var        0.005064   0.373402   0.014    0.989
beta0                -0.237692   0.460689  -0.516    0.606
beta5                -0.523521   0.456587  -1.147    0.252
cat_pre_wt_log_scale  0.330169   0.220049   1.500    0.134
car::Anova(eclosure_m, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: eclosed
                      Chisq Df Pr(>Chisq)
var_trt              0.0002  1     0.9892
beta                 1.3212  2     0.5165
cat_pre_wt_log_scale 2.2513  1     0.1335

2.3.2 Post-hoc contrasts

2.3.3 Misc stats

r2(eclosure_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.092
     Marginal R2: 0.046

2.3.4 Check model

check_mod(eclosure_m)

2.3.5 Plot

2.4 Survival

2.4.1 Model summary

surv_m_full <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m_full)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) * 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                                       coef exp(coef) se(coef)      z Pr(>|z|)
var_trtlow_var                      -0.1115    0.8945   0.2124 -0.525 0.599727
beta0                               -0.0896    0.9143   0.2579 -0.347 0.728306
beta5                               -0.2521    0.7772   0.2605 -0.968 0.333144
cat_pre_wt_log_scale                 0.9641    2.6224   0.2556  3.772 0.000162
I(cat_pre_wt_log_scale^2)            0.2450    1.2776   0.1447  1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626    0.8499   0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale          -0.2230    0.8001   0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale          -0.4789    0.6195   0.2584 -1.854 0.063807
                                       
var_trtlow_var                         
beta0                                  
beta5                                  
cat_pre_wt_log_scale                ***
I(cat_pre_wt_log_scale^2)           .  
var_trtlow_var:cat_pre_wt_log_scale    
beta0:cat_pre_wt_log_scale             
beta5:cat_pre_wt_log_scale          .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                                    exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var                         0.8945     1.1179    0.5899     1.356
beta0                                  0.9143     1.0937    0.5515     1.516
beta5                                  0.7772     1.2867    0.4665     1.295
cat_pre_wt_log_scale                   2.6224     0.3813    1.5892     4.328
I(cat_pre_wt_log_scale^2)              1.2776     0.7827    0.9621     1.696
var_trtlow_var:cat_pre_wt_log_scale    0.8499     1.1766    0.5746     1.257
beta0:cat_pre_wt_log_scale             0.8001     1.2499    0.5089     1.258
beta5:cat_pre_wt_log_scale             0.6195     1.6142    0.3734     1.028

Concordance= 0.637  (se = 0.047 )
Likelihood ratio test= 22  on 8 df,   p=0.005
Wald test            = 22.31  on 8 df,   p=0.004
Score (logrank) test = 24.57  on 8 df,   p=0.002
drop1(surv_m_full, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                             Df    AIC    LRT Pr(>Chi)  
<none>                          459.50                  
I(cat_pre_wt_log_scale^2)     1 460.40   2.90  0.08841 .
strata(session_id)            0 787.87 328.37           
var_trt:cat_pre_wt_log_scale  1 458.16   0.67  0.41395  
beta:cat_pre_wt_log_scale     2 458.98   3.48  0.17532  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
surv_m <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + 
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) + 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                              coef exp(coef) se(coef)      z Pr(>|z|)    
var_trtlow_var            -0.07077   0.93167  0.20845 -0.340 0.734212    
beta0                     -0.07410   0.92858  0.25808 -0.287 0.774026    
beta5                     -0.17722   0.83759  0.25804 -0.687 0.492206    
cat_pre_wt_log_scale       0.65947   1.93376  0.18164  3.631 0.000283 ***
I(cat_pre_wt_log_scale^2)  0.29584   1.34425  0.14124  2.095 0.036205 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                          exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var               0.9317     1.0733    0.6192     1.402
beta0                        0.9286     1.0769    0.5599     1.540
beta5                        0.8376     1.1939    0.5051     1.389
cat_pre_wt_log_scale         1.9338     0.5171    1.3545     2.761
I(cat_pre_wt_log_scale^2)    1.3443     0.7439    1.0192     1.773

Concordance= 0.602  (se = 0.047 )
Likelihood ratio test= 18.08  on 5 df,   p=0.003
Wald test            = 18.29  on 5 df,   p=0.003
Score (logrank) test = 19.42  on 5 df,   p=0.002
drop1(surv_m, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                          Df    AIC    LRT  Pr(>Chi)    
<none>                       457.42                     
var_trt                    1 455.53   0.12 0.7341150    
beta                       2 453.89   0.48 0.7883535    
cat_pre_wt_log_scale       1 469.05  13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2)  1 459.84   4.42 0.0355376 *  
strata(session_id)         0 784.70 327.28              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.4.2 Post-hoc contrasts

2.4.3 Misc stats

r2(surv_m, tolerance = 10^-10)
  Nagelkerke's R2: 0.135

2.4.4 Check model

cox.zph(surv_m)
                          chisq df    p
var_trt                    1.60  1 0.21
beta                       1.43  2 0.49
cat_pre_wt_log_scale       2.61  1 0.11
I(cat_pre_wt_log_scale^2)  1.05  1 0.31
GLOBAL                     5.53  5 0.35

2.4.5 Plot

3 Structural Equation Model

3.1 Sub-models

3.1.1 RGR

subm_rgr <- glmmTMB(
  RGR_scale ~ 
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq +
    mean_toxic_conc_scale + 
    area_herb_log_scale * cat_pre_wt_log_scale + 
    sl_mean_obs_log_scale + 
    var_toxic_12_scale +
    (1|session_id),
  data = d2
); summary(subm_rgr)
 Family: gaussian  ( identity )
Formula:          RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +  
    mean_toxic_conc_scale + area_herb_log_scale * cat_pre_wt_log_scale +  
    sl_mean_obs_log_scale + var_toxic_12_scale + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   140.2    165.6    -60.1    120.2       84 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01202  0.1096  
 Residual               0.20224  0.4497  
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.202 

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)
(Intercept)                               0.041726   0.078479   0.532  0.59494
cat_pre_wt_log_scale                     -3.367086   0.405069  -8.312  < 2e-16
cat_pre_wt_log_scale_sq                  -2.789806   0.417379  -6.684 2.32e-11
mean_toxic_conc_scale                    -0.233908   0.068928  -3.393  0.00069
area_herb_log_scale                       0.773436   0.091519   8.451  < 2e-16
sl_mean_obs_log_scale                    -0.164921   0.053712  -3.070  0.00214
var_toxic_12_scale                       -0.005495   0.055009  -0.100  0.92043
cat_pre_wt_log_scale:area_herb_log_scale -0.200719   0.091201  -2.201  0.02775
                                            
(Intercept)                                 
cat_pre_wt_log_scale                     ***
cat_pre_wt_log_scale_sq                  ***
mean_toxic_conc_scale                    ***
area_herb_log_scale                      ***
sl_mean_obs_log_scale                    ** 
var_toxic_12_scale                          
cat_pre_wt_log_scale:area_herb_log_scale *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.2 Variance in toxin (12 hour)

subm_var_toxic <- glmmTMB(
  var_toxic_12_scale ~ 
    beta_numeric_scale + 
    var_high + 
    sl_mean_obs_log_scale + 
    select_scale + 
    ava_qual_logit_scale + 
    (1|session_id),
  data = d2
); summary(subm_var_toxic)
 Family: gaussian  ( identity )
Formula:          
var_toxic_12_scale ~ beta_numeric_scale + var_high + sl_mean_obs_log_scale +  
    select_scale + ava_qual_logit_scale + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   163.9    184.2    -73.9    147.9       86 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01202  0.1097  
 Residual               0.27355  0.5230  
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.274 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -0.86243    0.09523  -9.057  < 2e-16 ***
beta_numeric_scale    -0.10881    0.05830  -1.866 0.062004 .  
var_high               1.69712    0.11500  14.757  < 2e-16 ***
sl_mean_obs_log_scale  0.13145    0.06474   2.030 0.042331 *  
select_scale          -0.20633    0.08159  -2.529 0.011440 *  
ava_qual_logit_scale  -0.28324    0.07880  -3.594 0.000325 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.3 Mean step length

subm_sl <- glmmTMB(
  sl_mean_obs_log_scale ~ 
    beta_numeric_scale + 
    on_toxic_logit_scale + 
    var_high +  
    select_scale + 
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq + 
    (1|session_id),
  data = d2,
); summary(subm_sl)
 Family: gaussian  ( identity )
Formula:          
sl_mean_obs_log_scale ~ beta_numeric_scale + on_toxic_logit_scale +  
    var_high + select_scale + cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   245.0    267.9   -113.5    227.0       85 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.828e-10 1.682e-05
 Residual               6.551e-01 8.094e-01
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.655 

Conditional model:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)              0.24784    0.12625   1.963 0.049627 *  
beta_numeric_scale       0.18760    0.08691   2.159 0.030879 *  
on_toxic_logit_scale     0.38465    0.11614   3.312 0.000926 ***
var_high                -0.24283    0.17679  -1.374 0.169587    
select_scale             0.35697    0.13067   2.732 0.006297 ** 
cat_pre_wt_log_scale     1.51783    0.58084   2.613 0.008971 ** 
cat_pre_wt_log_scale_sq  1.84748    0.61197   3.019 0.002537 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.4 Mean toxin ingested

subm_toxin_ingested <- glmmTMB(
  mean_toxic_conc_scale ~ 
    beta_numeric_scale + 
    var_high * cat_pre_wt_log_scale + 
    on_toxic_logit_scale + 
    cat_pre_wt_log_scale_sq +  
    (1|session_id),
  data = d2 
); summary(subm_toxin_ingested)
 Family: gaussian  ( identity )
Formula:          
mean_toxic_conc_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +  
    on_toxic_logit_scale + cat_pre_wt_log_scale_sq + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   154.8    177.7    -68.4    136.8       85 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev. 
 session_id (Intercept) 6.90e-11 8.307e-06
 Residual               2.51e-01 5.010e-01
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.251 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.26244    0.07904   3.320 0.000899 ***
beta_numeric_scale            -0.08727    0.05331  -1.637 0.101636    
var_high                      -0.81255    0.11198  -7.257 3.97e-13 ***
cat_pre_wt_log_scale           0.96560    0.37021   2.608 0.009100 ** 
on_toxic_logit_scale           0.23369    0.06094   3.835 0.000126 ***
cat_pre_wt_log_scale_sq        1.06006    0.37862   2.800 0.005113 ** 
var_high:cat_pre_wt_log_scale -0.26558    0.11647  -2.280 0.022596 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.5 Mean time on toxin

subm_on_toxic <- glmmTMB(
  on_toxic_logit_scale ~ 
    select_scale + ava_qual_logit_scale  + 
    (1|session_id),
  family = gaussian(),
  data = d2
); summary(subm_on_toxic)
 Family: gaussian  ( identity )
Formula:          on_toxic_logit_scale ~ select_scale + ava_qual_logit_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
  -106.0    -93.3     58.0   -116.0       89 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 9.663e-12 3.109e-06
 Residual               1.705e-02 1.306e-01
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.0171 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -0.009281   0.013505   -0.69    0.492    
select_scale         -0.157732   0.019369   -8.14 3.84e-16 ***
ava_qual_logit_scale -0.909758   0.016169  -56.27  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.6 Area consumption

subm_herb <- glmmTMB(
  area_herb_log_scale ~ 
    beta_numeric_scale + var_high * cat_pre_wt_log_scale + 
    (1|session_id),
  data = d2,
); summary(subm_herb)
 Family: gaussian  ( identity )
Formula:          
area_herb_log_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   146.4    164.2    -66.2    132.4       87 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.742e-11 8.211e-06
 Residual               2.394e-01 4.893e-01
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.239 

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.046979   0.074660   0.629  0.52919    
beta_numeric_scale             0.001622   0.050860   0.032  0.97457    
var_high                       0.283395   0.104609   2.709  0.00675 ** 
cat_pre_wt_log_scale           0.550474   0.087850   6.266  3.7e-10 ***
var_high:cat_pre_wt_log_scale -0.355786   0.113217  -3.143  0.00168 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.7 Toxin selection

subm_select <- glmmTMB(
  select_scale ~ 
    beta_numeric_scale + var_high * cat_pre_wt_log_scale + 
    ava_qual_logit_scale + 
    (1|session_id),
  data = d2,
); summary(subm_select)
 Family: gaussian  ( identity )
Formula:          
select_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +  
    ava_qual_logit_scale + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   205.9    226.3    -95.0    189.9       86 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 1.243e-10 1.115e-05
 Residual               4.415e-01 6.645e-01
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.442 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -0.03222    0.10287  -0.313   0.7541    
beta_numeric_scale            -0.08040    0.07104  -1.132   0.2577    
var_high                      -0.01374    0.14649  -0.094   0.9253    
cat_pre_wt_log_scale          -0.12953    0.12030  -1.077   0.2816    
ava_qual_logit_scale           0.37617    0.07942   4.736 2.18e-06 ***
var_high:cat_pre_wt_log_scale  0.37178    0.15379   2.417   0.0156 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.8 Neighborhood diet quality

subm_ava <- glmmTMB(
  ava_qual_logit_scale ~ 
    var_high + beta_numeric_scale * cat_pre_wt_log_scale + 
    (1|session_id),
  family = gaussian(),
  data = d2,
); summary(subm_ava)
 Family: gaussian  ( identity )
Formula:          
ava_qual_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   248.3    266.1   -117.1    234.3       87 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01866  0.1366  
 Residual               0.69278  0.8323  
Number of obs: 94, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.693 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)   
(Intercept)                             -0.20763    0.14063  -1.476  0.13983   
var_high                                 0.45433    0.17311   2.625  0.00868 **
beta_numeric_scale                       0.15703    0.09021   1.741  0.08173 . 
cat_pre_wt_log_scale                     0.18432    0.12145   1.518  0.12912   
beta_numeric_scale:cat_pre_wt_log_scale  0.21360    0.09834   2.172  0.02985 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2 d-separation test

sem_fit <- psem(
    subm_rgr,
    subm_var_toxic,
    subm_on_toxic,
    subm_toxin_ingested,
    subm_select,
    subm_ava,
    subm_herb,
    subm_sl, 
    var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
    data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
  summary_psem(sem_fit, 
               no_standardize_x =  c("var_high","beta_red","beta_white"), 
               .progressBar = FALSE, 
               # Remove 'ava_qual_logit_scale' from the independence claim when 'on_toxic_logit_scale' is conditioned. 'on_toxic_logit_scale' is a more proximal predictor and is highly collinear with ava_qual_logit_scale. 
               basis.set = keep(basisSet(sem_fit), function(x){ 
                 res1 <- (basisSet_predictor(x) == "ava_qual_logit_scale")
                 res2 <- ("on_toxic_logit_scale" %in% basisSet_conditioned(x))
                 res <- res1 & res2
                 return(!res)
               }))
  ))
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
sem_summary$Cstat
  Fisher.C df P.Value
1   48.845 48   0.439
knitr::kable(sem_summary$dTable)
Independ.Claim Test.Type DF Crit.Value P.Value
var_toxic_12_scale ~ cat_pre_wt_log_scale + … coef 94 1.5236 0.1276
on_toxic_logit_scale ~ cat_pre_wt_log_scale + … coef 94 0.1994 0.8419
area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … coef 94 0.8087 0.4187
var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … coef 94 -1.2587 0.2081
select_scale ~ cat_pre_wt_log_scale_sq + … coef 94 -0.8184 0.4131
ava_qual_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 94 -0.4761 0.6340
on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 94 -0.3624 0.7171
RGR_scale ~ beta_numeric_scale + … coef 94 -1.8383 0.0660
on_toxic_logit_scale ~ beta_numeric_scale + … coef 94 0.1477 0.8826
RGR_scale ~ var_high + … coef 94 -0.3104 0.7563
on_toxic_logit_scale ~ var_high + … coef 94 -0.6569 0.5113
select_scale ~ RGR_scale + … coef 94 -1.0670 0.2860
ava_qual_logit_scale ~ RGR_scale + … coef 94 1.6016 0.1093
on_toxic_logit_scale ~ RGR_scale + … coef 94 -0.7917 0.4285
sl_mean_obs_log_scale ~ area_herb_log_scale + … coef 94 -0.3715 0.7103
var_toxic_12_scale ~ area_herb_log_scale + … coef 94 -1.5513 0.1208
select_scale ~ area_herb_log_scale + … coef 94 -0.0253 0.9798
ava_qual_logit_scale ~ area_herb_log_scale + … coef 94 1.0211 0.3072
on_toxic_logit_scale ~ area_herb_log_scale + … coef 94 1.1450 0.2522
mean_toxic_conc_scale ~ area_herb_log_scale + … coef 94 0.5116 0.6089
mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … coef 94 -0.4908 0.6236
mean_toxic_conc_scale ~ var_toxic_12_scale + … coef 94 -1.4198 0.1557
mean_toxic_conc_scale ~ select_scale + … coef 94 1.1730 0.2408
mean_toxic_conc_scale ~ ava_qual_logit_scale + … coef 94 0.6698 0.5030

3.3 Model summary

3.3.1 Model stats

sem_summary$ChiSq
NULL
sem_summary$AIC
       AIC     AICc  K  n
1 1198.514 1212.125 63 94
sem_summary$R2
               Response   family     link method Marginal Conditional
1             RGR_scale gaussian identity   none     0.74        0.76
2    var_toxic_12_scale gaussian identity   none     0.72        0.73
3  on_toxic_logit_scale gaussian identity   none     0.98          NA
4 mean_toxic_conc_scale gaussian identity   none     0.60          NA
5          select_scale gaussian identity   none     0.29          NA
6  ava_qual_logit_scale gaussian identity   none     0.19        0.21
7   area_herb_log_scale gaussian identity   none     0.36          NA
8 sl_mean_obs_log_scale gaussian identity   none     0.26          NA

3.3.2 Standardized coefficients

knitr::kable(
  cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*10))
)
Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate star arrow_size
RGR_scale cat_pre_wt_log_scale -3.3671 0.4051 94 -8.3124 0.0000 -3.4487000 *** 34.4870000
RGR_scale cat_pre_wt_log_scale_sq -2.7898 0.4174 94 -6.6841 0.0000 -2.7190000 *** 27.1900000
RGR_scale mean_toxic_conc_scale -0.2339 0.0689 94 -3.3935 0.0007 -0.2088000 *** 2.0880000
RGR_scale area_herb_log_scale 0.7734 0.0915 94 8.4511 0.0000 0.5289000 *** 5.2890000
RGR_scale sl_mean_obs_log_scale -0.1649 0.0537 94 -3.0704 0.0021 -0.1739000 ** 1.7390000
RGR_scale var_toxic_12_scale -0.0055 0.055 94 -0.0999 0.9204 -0.0061000 0.0610000
RGR_scale cat_pre_wt_log_scale:area_herb_log_scale -0.2007 0.0912 94 -2.2008 0.0277 -0.1353000 * 1.3530000
var_toxic_12_scale beta_numeric_scale -0.1088 0.0583 94 -1.8663 0.0620 -0.1095000 1.0950000
var_toxic_12_scale var_high 1.6971 0.115 94 14.7574 0.0000 1.7241172 *** 17.2411719
var_toxic_12_scale sl_mean_obs_log_scale 0.1314 0.0647 94 2.0303 0.0423 0.1248000 * 1.2480000
var_toxic_12_scale select_scale -0.2063 0.0816 94 -2.5290 0.0114 -0.1643000 * 1.6430000
var_toxic_12_scale ava_qual_logit_scale -0.2832 0.0788 94 -3.5942 0.0003 -0.2702000 *** 2.7020000
on_toxic_logit_scale select_scale -0.1577 0.0194 94 -8.1434 0.0000 -0.1333000 *** 1.3330000
on_toxic_logit_scale ava_qual_logit_scale -0.9098 0.0162 94 -56.2663 0.0000 -0.9207000 *** 9.2070000
mean_toxic_conc_scale beta_numeric_scale -0.0873 0.0533 94 -1.6370 0.1016 -0.1093000 1.0930000
mean_toxic_conc_scale var_high -0.8126 0.112 94 -7.2565 0.0000 -1.0180231 *** 10.1802310
mean_toxic_conc_scale cat_pre_wt_log_scale 0.9656 0.3702 94 2.6083 0.0091 1.1079000 ** 11.0790000
mean_toxic_conc_scale on_toxic_logit_scale 0.2337 0.0609 94 3.8346 0.0001 0.2741000 *** 2.7410000
mean_toxic_conc_scale cat_pre_wt_log_scale_sq 1.0601 0.3786 94 2.7998 0.0051 1.1574000 ** 11.5740000
mean_toxic_conc_scale var_high:cat_pre_wt_log_scale -0.2656 0.1165 94 -2.2802 0.0226 -0.2404000 * 2.4040000
select_scale beta_numeric_scale -0.0804 0.071 94 -1.1318 0.2577 -0.1016000 1.0160000
select_scale var_high -0.0137 0.1465 94 -0.0938 0.9253 -0.0173199 0.1731994
select_scale cat_pre_wt_log_scale -0.1295 0.1203 94 -1.0768 0.2816 -0.1500000 1.5000000
select_scale ava_qual_logit_scale 0.3762 0.0794 94 4.7364 0.0000 0.4506000 *** 4.5060000
select_scale var_high:cat_pre_wt_log_scale 0.3718 0.1538 94 2.4174 0.0156 0.3396000 * 3.3960000
ava_qual_logit_scale var_high 0.4543 0.1731 94 2.6246 0.0087 0.5111900 ** 5.1119004
ava_qual_logit_scale beta_numeric_scale 0.1570 0.0902 94 1.7407 0.0817 0.1657000 1.6570000
ava_qual_logit_scale cat_pre_wt_log_scale 0.1843 0.1215 94 1.5176 0.1291 0.1782000 1.7820000
ava_qual_logit_scale beta_numeric_scale:cat_pre_wt_log_scale 0.2136 0.0983 94 2.1721 0.0298 0.2069000 * 2.0690000
area_herb_log_scale beta_numeric_scale 0.0016 0.0509 94 0.0319 0.9746 0.0027000 0.0270000
area_herb_log_scale var_high 0.2834 0.1046 94 2.7091 0.0067 0.4634544 ** 4.6345440
area_herb_log_scale cat_pre_wt_log_scale 0.5505 0.0878 94 6.2661 0.0000 0.8245000 *** 8.2450000
area_herb_log_scale var_high:cat_pre_wt_log_scale -0.3558 0.1132 94 -3.1425 0.0017 -0.4204000 ** 4.2040000
sl_mean_obs_log_scale beta_numeric_scale 0.1876 0.0869 94 2.1586 0.0309 0.1989000 * 1.9890000
sl_mean_obs_log_scale on_toxic_logit_scale 0.3847 0.1161 94 3.3120 0.0009 0.3819000 *** 3.8190000
sl_mean_obs_log_scale var_high -0.2428 0.1768 94 -1.3735 0.1696 -0.2574872 2.5748724
sl_mean_obs_log_scale select_scale 0.3570 0.1307 94 2.7319 0.0063 0.2994000 ** 2.9940000
sl_mean_obs_log_scale cat_pre_wt_log_scale 1.5178 0.5808 94 2.6131 0.0090 1.4742000 ** 14.7420000
sl_mean_obs_log_scale cat_pre_wt_log_scale_sq 1.8475 0.612 94 3.0189 0.0025 1.7075000 ** 17.0750000
~~var_toxic_12_scale ~~on_toxic_logit_scale 0.0823 - 94 0.7882 0.2163 0.0823000 0.8230000

3.3.3 Post hoc contrast

contrast_by_pre_wt(subm_toxin_ingested, "var_high", "trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_high emmean    SE df lower.CL upper.CL
        0 -1.941 0.823 85   -3.577   -0.306
        1 -2.223 0.819 85   -3.851   -0.594

cat_pre_wt_log_scale =  2:
 var_high emmean    SE df lower.CL upper.CL
        0  1.921 0.666 85    0.597    3.246
        1  0.578 0.629 85   -0.672    1.827

Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -0.281 0.279 85  -1.008  0.3163

cat_pre_wt_log_scale =  2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -1.344 0.236 85  -5.695  <.0001
1/sd(subm_toxin_ingested$frame$mean_toxic_conc_scale)
[1] 1.252797
contrast_by_pre_wt(subm_select, "var_high", "trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_high emmean    SE df lower.CL upper.CL
        0  0.251 0.287 86  -0.3197   0.8211
        1 -0.507 0.241 86  -0.9854  -0.0278

cat_pre_wt_log_scale =  2:
 var_high emmean    SE df lower.CL upper.CL
        0 -0.267 0.235 86  -0.7336   0.1988
        1  0.462 0.201 86   0.0638   0.8610

Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -0.757 0.371 86  -2.041  0.0443

cat_pre_wt_log_scale =  2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0    0.730 0.307 86   2.375  0.0198
1/sd(subm_select$frame$select_scale)
[1] 1.264229
contrast_by_pre_wt(subm_herb, "var_high", "trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_high emmean    SE df lower.CL upper.CL
        0 -1.054 0.208 87   -1.467   -0.641
        1 -0.059 0.176 87   -0.410    0.292

cat_pre_wt_log_scale =  2:
 var_high emmean    SE df lower.CL upper.CL
        0  1.148 0.173 87    0.805    1.491
        1  0.720 0.143 87    0.436    1.004

Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0    0.995 0.272 87   3.652  0.0004

cat_pre_wt_log_scale =  2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -0.428 0.224 87  -1.911  0.0593
1/sd(subm_herb$frame$area_herb_log_scale)
[1] 1.635337
emmeans::emtrends(subm_rgr, 
                 ~ area_herb_log_scale | cat_pre_wt_log_scale, 
                 var = "area_herb_log_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 area_herb_log_scale area_herb_log_scale.trend    SE df lower.CL upper.CL
               0.276                     1.175 0.190 84   0.7974    1.552

cat_pre_wt_log_scale =  2:
 area_herb_log_scale area_herb_log_scale.trend    SE df lower.CL upper.CL
               0.276                     0.372 0.217 84  -0.0603    0.804

Confidence level used: 0.95 
sd(subm_rgr$frame$area_herb_log_scale) / sd(subm_rgr$frame$RGR_scale)
[1] 0.6838414
emmeans::emtrends(subm_ava, 
                 ~ beta_numeric_scale | cat_pre_wt_log_scale, 
                 var = "beta_numeric_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.95e-17                   -0.270 0.238 87   -0.744    0.203

cat_pre_wt_log_scale =  2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.95e-17                    0.584 0.192 87    0.202    0.966

Results are averaged over the levels of: var_high 
Confidence level used: 0.95 
sd(subm_ava$frame$beta_numeric_scale) * sd(subm_ava$frame$ava_qual_logit_scale)
[1] 0.9475721

3.4 SEM plots

knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_joined.png"), rel_path = FALSE)

3.5 SEM forest plot

knitr::include_graphics(paste0(getwd(),"/graphs/manuscript1_figures/SEM_forest.png"), rel_path = FALSE)

3.6 Misc

3.6.1 Time to pupation

subm_pupt <- glmmTMB(
  time_to_pupation ~ 
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq +
    mean_toxic_conc_scale + 
    area_herb_log_scale  + 
    sl_mean_obs_log_scale + 
    var_toxic_12_scale +
    (1|session_id),
  family = Gamma(link = "log"),
  data = d2 %>% 
    filter(!is.na(time_to_pupation))
); summary(subm_pupt)
 Family: Gamma  ( log )
Formula:          
time_to_pupation ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +  
    mean_toxic_conc_scale + area_herb_log_scale + sl_mean_obs_log_scale +  
    var_toxic_12_scale + (1 | session_id)
Data: d2 %>% filter(!is.na(time_to_pupation))

     AIC      BIC   logLik deviance df.resid 
   280.8    302.8   -131.4    262.8       76 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 3.69e-05 0.006075
Number of obs: 85, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0224 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)              2.205884   0.021976  100.38  < 2e-16 ***
cat_pre_wt_log_scale    -1.107727   0.137192   -8.07 6.79e-16 ***
cat_pre_wt_log_scale_sq -0.781595   0.150144   -5.21 1.93e-07 ***
mean_toxic_conc_scale    0.082351   0.024635    3.34 0.000829 ***
area_herb_log_scale     -0.139853   0.044549   -3.14 0.001693 ** 
sl_mean_obs_log_scale    0.058968   0.019041    3.10 0.001956 ** 
var_toxic_12_scale       0.003529   0.020361    0.17 0.862383    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(subm_pupt)
# R2 for Mixed Models

  Conditional R2: 0.888
     Marginal R2: 0.888

4 Session Info

sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] spat1f_0.0.1       herbivar_0.2.1     forcats_0.5.1      stringr_1.5.0     
 [5] dplyr_1.1.2        purrr_1.0.1        readr_2.1.2        tidyr_1.3.0       
 [9] tibble_3.2.1       tidyverse_1.3.1    imager_0.42.13     magrittr_2.0.3    
[13] doSNOW_1.0.20      snow_0.4-4         iterators_1.0.14   foreach_1.5.2     
[17] glmmTMB_1.1.7      tidybayes_3.0.2    performance_0.10.4 sjPlot_2.8.15     
[21] ggpubr_0.6.0       ggplot2_3.4.2      piecewiseSEM_2.3.0 survival_3.3-1    

loaded via a namespace (and not attached):
  [1] svUnit_1.0.6             splines_4.3.1            later_1.3.0             
  [4] cellranger_1.1.0         polyclip_1.10-0          reprex_2.0.1            
  [7] lifecycle_1.0.3          sf_1.0-8                 Rdpack_2.3.1            
 [10] rstatix_0.7.2            doParallel_1.0.17        rprojroot_2.0.3         
 [13] vroom_1.5.7              processx_3.6.1           lattice_0.21-8          
 [16] MASS_7.3-60              insight_0.19.3           readbitmap_0.1.5        
 [19] ggdist_3.1.1             backports_1.4.1          sass_0.4.1              
 [22] rmarkdown_2.24           jquerylib_0.1.4          yaml_2.3.5              
 [25] remotes_2.4.2            httpuv_1.6.5             gap_1.2.3-6             
 [28] sessioninfo_1.2.2        pkgbuild_1.3.1           cowplot_1.1.1           
 [31] DBI_1.1.3                minqa_1.2.4              RColorBrewer_1.1-3      
 [34] lubridate_1.8.0          DHARMa_0.4.6             multcomp_1.4-25         
 [37] abind_1.4-5              pkgload_1.3.2            rvest_1.0.2             
 [40] imagerExtra_2.0.0        TH.data_1.1-1            tensorA_0.36.2          
 [43] sandwich_3.0-2           spatstat.utils_3.0-1     units_0.8-0             
 [46] codetools_0.2-19         MuMIn_1.47.5             xml2_1.3.3              
 [49] tidyselect_1.2.0         ggeffects_1.2.3          farver_2.1.0            
 [52] lme4_1.1-29              matrixStats_1.2.0        stats4_4.3.1            
 [55] jsonlite_1.8.8           e1071_1.7-11             ellipsis_0.3.2          
 [58] bmp_0.3                  emmeans_1.7.5            tools_4.3.1             
 [61] Rcpp_1.0.8.3             glue_1.6.2               mgcv_1.8-40             
 [64] xfun_0.39                distributional_0.3.0     usethis_2.1.6           
 [67] withr_2.5.0              numDeriv_2016.8-1.1      fastmap_1.1.0           
 [70] boot_1.3-28.1            fansi_1.0.3              callr_3.7.0             
 [73] digest_0.6.29            R6_2.5.1                 mime_0.12               
 [76] estimability_1.3         colorspace_2.0-3         spatstat.data_3.0-0     
 [79] jpeg_0.1-9               see_0.7.1                DiagrammeR_1.0.9        
 [82] utf8_1.2.2               generics_0.1.2           class_7.3-22            
 [85] prettyunits_1.1.1        httr_1.4.3               htmlwidgets_1.5.4       
 [88] pkgconfig_2.0.3          gtable_0.3.0             htmltools_0.5.2         
 [91] carData_3.0-5            profvis_0.3.7            fftwtools_0.9-11        
 [94] TMB_1.9.4                scales_1.2.1             png_0.1-7               
 [97] posterior_1.2.2          knitr_1.43               rstudioapi_0.13         
[100] tzdb_0.3.0               coda_0.19-4              visNetwork_2.1.0        
[103] checkmate_2.1.0          nlme_3.1-162             nloptr_2.0.3            
[106] proxy_0.4-27             cachem_1.0.6             zoo_1.8-12              
[109] sjlabelled_1.2.0         KernSmooth_2.23-20       parallel_4.3.1          
[112] miniUI_0.1.1.1           desc_1.4.1               pillar_1.9.0            
[115] grid_4.3.1               vctrs_0.6.3              urlchecker_1.0.1        
[118] promises_1.2.0.1         car_3.1-2                dbplyr_2.2.0            
[121] arrayhelpers_1.1-0       xtable_1.8-4             evaluate_0.15           
[124] amt_0.2.1.0              mvtnorm_1.1-3            cli_3.6.1               
[127] compiler_4.3.1           rlang_1.1.1              crayon_1.5.1            
[130] ggsignif_0.6.3           labeling_0.4.2           modelr_0.1.8            
[133] classInt_0.4-7           RcppArmadillo_0.11.2.0.0 ps_1.7.1                
[136] plyr_1.8.7               sjmisc_2.8.9             fs_1.5.2                
[139] gap.datasets_0.0.5       stringi_1.7.12           deldir_1.0-6            
[142] assertthat_0.2.1         munsell_0.5.0            tiff_0.1-11             
[145] spatstat.geom_3.0-3      devtools_2.4.5           bayestestR_0.13.1       
[148] Matrix_1.5-4.1           sjstats_0.18.2           hms_1.1.1               
[151] bit64_4.0.5              qgam_1.3.4               shiny_1.7.1             
[154] highr_0.9                haven_2.5.0              rbibutils_2.2.8         
[157] igraph_1.3.2             broom_1.0.5              memoise_2.0.1           
[160] bslib_0.3.1              bit_4.0.4                readxl_1.4.0